home *** CD-ROM | disk | FTP | other *** search
/ McDisk 1 / McDisk 1.adf / sources / Floatpoint_root.S next >
Text File  |  1978-04-12  |  2KB  |  111 lines
  1. ;*************************************************
  2. ;* ROOT-ROUTINE FOR 32-BIT-FLOATINGPOINT-NUMBERS *
  3. ;*          coded by Otto Peter, Inbruck         *
  4. ;*             published in  C`T 1/90            *
  5. ;*************************************************
  6.  
  7. exp_offset:    equ $7f
  8.  
  9.  
  10. domain_error:            ;Error if number is negative
  11. movem.l    (a7)+,d1-d4
  12. rts
  13.  
  14. sq_0:
  15. clr.l    d0
  16. movem.l    (a7)+,d1-d4
  17. rts
  18.  
  19. sqrt:
  20. movem.l    d1-d4,-(a7)        ;Save d1-d4
  21. move.l    d0,d4
  22. bmi.s    domain_error        ;Error if number is negative
  23. swap    d4            ;MSW of number
  24. and.l    #$7f80,d0        ;isolate exponent
  25. beq.s    sq_0            ;exponent is 0 => root is 0
  26. and.l    #$007fffff,d0        ;isolate mantisse
  27. sub.w    #exp_offset*$80,d4    ;exponent in `2^i`-format
  28. bclr    #7,d4            ;exponent even?
  29. beq.s    even_exp
  30. add.l    d0,d0            ;mantisse * 2
  31. add.l    #$01000000-$00800000,d0    ;set hidden-bit
  32.  
  33. even_exp:
  34. asr.w    #1,d4            ;exponent / 2
  35. add.w    #exp_offset*$80,d4    ;exponent in offset-format
  36. swap    d4
  37.  
  38. lsl.l    #7,d0
  39.  
  40. move.l    #$40000000,d2        ;xroot after 1. interation
  41. move.l    #$10000000,d3        ;m2 = 2 << (maxbit-1)
  42.  
  43. loop1_0:
  44. move.l    d0,d1            ;xx2 = x
  45.  
  46. loop1_1:
  47. sub.l    d2,d1            ;xx2 -= xroot
  48. lsr.l    #1,d2            ;xroot >>= 1
  49. sub.l    d3,d1            ;x2 -= m2
  50. bmi.s    dont_set1
  51. move.l    d1,d0            ;x = xx2
  52. or.l    d3,d2            ;xroot += m2
  53. lsr.l    #2,d3            ;m2 >>= 2
  54. bne.s    loop1_1
  55. bra.s    d0_d1_same
  56.  
  57. dont_set1:
  58. lsr.l    d2,d3            ;m2 >>= 2
  59. bne.s    loop1_0            ;repeat loop 15 times (bit 22..8)
  60.  
  61. move.l    d0,d1
  62.  
  63. d0_d1_same:
  64. sub.l    d2,d1            ;xx2 -= xroot
  65. ror.l    #1,d2            ;xroot >>= 1 (with carry)
  66. swap     d2            ;turn to new aligment
  67. subq.l    #1,d1            ;carry of 0-0x4000: x2 -= m2 (part 1)
  68. bmi.s    dont_set7
  69.  
  70. or.l    #-$40000000,d1        ;0-0x4000: x2 -= m2 (part 2)
  71. move.l    d1,d0
  72. or.w    #$4000,d2
  73.  
  74. dont_set7:
  75. swap    d0            ;turn x to new aligment
  76.  
  77. move.w    #$1000,d3        ;m2 - bit 16..31 = 0
  78.  
  79. loop2_0:
  80. move.l    d0,d1            ;xx2 = x
  81.  
  82. loop2_1:
  83. sub.l    d2,d1            ;xx2 = xroot
  84. lsr.l    #1,d2            ;xroot >>= 1
  85. sub.l    d3,d1
  86. bmi.s    dont_set2
  87.  
  88. move.l    d1,d0            ;x = xx2
  89. or.l    d3,d2            ;xroot += m2
  90. lsr.w    #2,d3            ;m2 >>= 2
  91. bne.s    loop2_1
  92.  
  93. bra.s    finish
  94.  
  95. dont_set2:
  96. lsr.w    #2,d3            ;m2 >>= 2
  97. bne.s    loop2_0            ;repeat loop 7 times (n=6..0)
  98.  
  99. finish:
  100. sub.l    d2,d0            ;round root?
  101. bls.s    no_inc
  102. addq.l    #1,d2            ;round up!
  103.  
  104. no_inc:
  105. bclr    #23,d2            ;clear hidden bit
  106. or.l    d4,d2            ;combine exponent and mantisse
  107. move.l    d2,d0            ;result
  108. movem.l    (a7)+,d1-d4
  109. rts
  110.  
  111.